module parallel_mod

	use const_mod
	use mesh_mod

	implicit none

	private

  public fill_zonal_halo_cell
	public fill_merid_halo_cell
  public fill_zonal_halo_lat
  public fill_zonal_halo_lon

  interface fill_zonal_halo_cell ! cell center
    module procedure fill_zonal_halo_cell_2d
    module procedure fill_zonal_halo_cell_3d
  end interface fill_zonal_halo_cell

  interface fill_merid_halo_cell ! cell center
    module procedure fill_merid_halo_cell_2d
    module procedure fill_merid_halo_cell_3d
  end interface

  interface fill_zonal_halo_lat ! lat
    module procedure fill_zonal_halo_lat_2d
    module procedure fill_zonal_halo_lat_3d
  end interface

  interface fill_zonal_halo_lon ! lon
    module procedure fill_zonal_halo_lon_2d
    module procedure fill_zonal_halo_lon_3d
  end interface

contains 

  subroutine fill_zonal_halo_cell_2d(field, all_halo, left_halo, right_halo)

  	real(r8), intent(inout) :: field(mesh%full_lon_lb:mesh%full_lon_ub,&
                                     mesh%full_lat_lb:mesh%full_lat_ub)
  	logical, intent(in), optional :: all_halo
  	logical, intent(in), optional :: left_halo
  	logical, intent(in), optional :: right_halo

  	integer i, j, m, n
    logical all_, left_, right_

    if (present(all_halo)) then
      all_ = all_halo
    else
      all_ = .true.
    end if
    if (present(left_halo)) then
      left_ = left_halo
    else
      left_ = .true.
    end if
    if (present(right_halo)) then
      right_ = right_halo
    else
      right_ = .true.
    end if

		if (all_ .or. left_) then
			! left halo
			m = lbound(field, 1) - 1
			n = ubound(field, 1) - 2 * mesh%lon_halo_width
			do j = mesh%full_lat_ibeg+1, mesh%full_lat_iend-1
				do i = 1, mesh%lon_halo_width
					field(m+i,j) = field(n+i,j)
				end do
			end do
		end if

		if (all_ .or. right_) then
		  ! right halo
		  m = ubound(field, 1) - mesh%lon_halo_width
		  n = lbound(field, 1) + mesh%lon_halo_width - 1
		  do j = mesh%full_lat_ibeg+1, mesh%full_lat_iend-1
		  	do i = 1, mesh%lon_halo_width
		  		field(m+i,j) = field(n+i,j)
		  	end do
		  end do
		end if 

  end subroutine fill_zonal_halo_cell_2d

  subroutine fill_zonal_halo_cell_3d(field, all_halo, left_halo, right_halo)

    real(r8), intent(inout) :: field(mesh%full_lon_lb:mesh%full_lon_ub, &
                                     mesh%full_lat_lb:mesh%full_lat_ub, &
                                     mesh%full_lev_lb:mesh%full_lev_ub)
    logical, intent(in), optional :: all_halo
    logical, intent(in), optional :: left_halo
    logical, intent(in), optional :: right_halo

    integer i, j, k, m, n
    logical all_, left_, right_

    if (present(all_halo)) then
      all_ = all_halo
    else
      all_ = .false.
    end if
    if (present(left_halo)) then
      left_ = left_halo
    else
      left_ = .false.
    end if
    if (present(right_halo)) then
      right_ = right_halo
    else
      right_ = .false.
    end if

    if (all_ .or. left_) then
      ! left halo
      m = lbound(field, 1) - 1
      n = ubound(field, 1) - 2 * mesh%lon_halo_width
      do k = mesh%full_lev_lb, mesh%full_lev_ub
        do j = mesh%full_lat_ibeg+1, mesh%full_lat_iend-1
          do i = 1, mesh%lon_halo_width
            field(m+i,j,k) = field(n+i,j,k)
          end do
        end do
      end do
    end if

    if (all_ .or. right_) then
      ! right halo
      m = ubound(field, 1) - mesh%lon_halo_width
      n = lbound(field, 1) + mesh%lon_halo_width - 1
      do k = mesh%full_lev_lb, mesh%full_lev_ub
        do j = mesh%full_lat_ibeg+1, mesh%full_lat_iend-1
          do i = 1, mesh%lon_halo_width
            field(m+i,j,k) = field(n+i,j,k)
          end do
        end do
      end do
    end if 

  end subroutine fill_zonal_halo_cell_3d

  subroutine fill_zonal_halo_lat_2d(field, all_halo, left_halo, right_halo)

  	real(r8), intent(inout) :: field(mesh%full_lon_lb:mesh%full_lon_ub, &
                                     mesh%half_lat_lb:mesh%half_lat_ub)
  	logical, intent(in), optional :: all_halo
  	logical, intent(in), optional :: left_halo
  	logical, intent(in), optional :: right_halo

  	integer i, j, m, n
    logical all_, left_, right_

    if (present(all_halo)) then
      all_ = all_halo
    else
      all_ = .true.
    end if
    if (present(left_halo)) then
      left_ = left_halo
    else
      left_ = .true.
    end if
    if (present(right_halo)) then
      right_ = right_halo
    else
      right_ = .true.
    end if
		if (all_ .or. left_) then
			! left halo
			m = lbound(field, 1) - 1
			n = ubound(field, 1) - 2 * mesh%lon_halo_width
			do j = mesh%half_lat_ibeg, mesh%half_lat_iend
				do i = 1, mesh%lon_halo_width
					field(m+i,j) = field(n+i,j)
				end do
			end do
		end if

		if (all_ .or. right_) then
		  ! right halo
		  m = ubound(field, 1) - mesh%lon_halo_width
		  n = lbound(field, 1) + mesh%lon_halo_width - 1
		  do j = mesh%half_lat_ibeg, mesh%half_lat_iend
		  	do i = 1, mesh%lon_halo_width
		  		field(m+i,j) = field(n+i,j)
		  	end do
		  end do
		end if 

  end subroutine fill_zonal_halo_lat_2d

  subroutine fill_zonal_halo_lat_3d(field, all_halo, left_halo, right_halo)

    real(r8), intent(inout) :: field(mesh%full_lon_lb:mesh%full_lon_ub, &
                                     mesh%half_lat_lb:mesh%half_lat_ub, &
                                     mesh%full_lev_lb:mesh%full_lev_ub)
    logical, intent(in), optional :: all_halo
    logical, intent(in), optional :: left_halo
    logical, intent(in), optional :: right_halo

    integer i, j, k, m, n
    logical all_, left_, right_

    if (present(all_halo)) then
      all_ = all_halo
    else
      all_ = .true.
    end if
    if (present(left_halo)) then
      left_ = left_halo
    else
      left_ = .true.
    end if
    if (present(right_halo)) then
      right_ = right_halo
    else
      right_ = .true.
    end if
    if (all_ .or. left_) then
      ! left halo
      m = lbound(field, 1) - 1
      n = ubound(field, 1) - 2 * mesh%lon_halo_width
      do k = mesh%full_lev_lb, mesh%full_lev_ub
        do j = mesh%half_lat_ibeg, mesh%half_lat_iend
          do i = 1, mesh%lon_halo_width
            field(m+i,j,k) = field(n+i,j,k)
          end do
        end do
      end do
    end if

    if (all_ .or. right_) then
      ! right halo
      m = ubound(field, 1) - mesh%lon_halo_width
      n = lbound(field, 1) + mesh%lon_halo_width - 1
      do k = mesh%full_lev_lb, mesh%full_lev_ub
        do j = mesh%half_lat_ibeg, mesh%half_lat_iend
          do i = 1, mesh%lon_halo_width
            field(m+i,j,k) = field(n+i,j,k)
          end do
        end do
      end do
    end if 

  end subroutine fill_zonal_halo_lat_3d

  subroutine fill_merid_halo_cell_2d(field, all_halo, south_halo, north_halo)

  	real(r8), intent(inout) :: field(mesh%full_lon_lb:mesh%full_lon_ub, &
                                     mesh%full_lat_lb:mesh%full_lat_ub)
  	logical, intent(in), optional :: all_halo
  	logical, intent(in), optional :: south_halo
  	logical, intent(in), optional :: north_halo

  	integer i, j, iop, m, n
    logical all_, south_, north_

    if (present(all_halo)) then
      all_ = all_halo
    else
      all_ = .true.
    end if

    if (present(south_halo)) then
      south_ = south_halo
    else
      south_ = .true.
    end if

    if (present(north_halo)) then
      north_ = north_halo
    else
      north_ = .true.
    end if
		if (all_ .or. south_) then
			n = lbound(field, 2) - 1
			do j = 1, mesh%lat_halo_width
  			do i = mesh%full_lon_ibeg, mesh%full_lon_iend
  				iop = i + mesh%num_full_lon / 2
  				if (iop > mesh%num_full_lon) then
  					iop = iop - mesh%num_full_lon 
  				end if    
  				field(i,n+j) = field(iop, 2 * mesh%lat_halo_width - j)
  			end do
  		end do
		end if


		if (all_ .or. north_) then
			n = ubound(field, 2) - mesh%lat_halo_width
			do j = 1, mesh%lat_halo_width
				do i = mesh%full_lon_ibeg, mesh%full_lon_iend
					iop = i + mesh%num_full_lon / 2
					if (iop > mesh%num_full_lon) then
						iop = iop - mesh%num_full_lon
					end if
					field(i,n+j) = field(iop, mesh%num_full_lat-j)
				end do
			end do
		end if

  end subroutine fill_merid_halo_cell_2d

  subroutine fill_merid_halo_cell_3d(field, all_halo, south_halo, north_halo)

    real(r8), intent(inout) :: field(mesh%full_lon_lb:mesh%full_lon_ub,&
                                     mesh%full_lat_lb:mesh%full_lat_ub,&
                                     mesh%full_lev_lb:mesh%full_lev_ub)
    logical, intent(in), optional :: all_halo
    logical, intent(in), optional :: south_halo
    logical, intent(in), optional :: north_halo

    integer i, j, k, iop, m, n
    logical all_, south_, north_

    if (present(all_halo)) then
      all_ = all_halo
    else
      all_ = .true.
    end if

    if (present(south_halo)) then
      south_ = south_halo
    else
      south_ = .true.
    end if

    if (present(north_halo)) then
      north_ = north_halo
    else
      north_ = .true.
    end if
    if (all_ .or. south_) then
      n = lbound(field, 2) - 1
      do k = mesh%full_lev_lb, mesh%full_lev_ub
        do j = 1, mesh%lat_halo_width
          do i = mesh%full_lon_ibeg, mesh%full_lon_iend
            iop = i + mesh%num_full_lon / 2
            if (iop > mesh%num_full_lon) then
              iop = iop - mesh%num_full_lon 
            end if    
            field(i,n+j,k) = field(iop,2 * mesh%lat_halo_width-j,k)
          end do
        end do
      end do
    end if


    if (all_ .or. north_) then
      n = ubound(field, 2) - mesh%lat_halo_width
      do k = mesh%full_lev_lb, mesh%full_lev_ub
        do j = 1, mesh%lat_halo_width
          do i = mesh%full_lon_ibeg, mesh%full_lon_iend
            iop = i + mesh%num_full_lon / 2
            if (iop > mesh%num_full_lon) then
              iop = iop - mesh%num_full_lon
            end if
            field(i,n+j,k) = field(iop,mesh%num_full_lat-j,k)
          end do
        end do
      end do
    end if

  end subroutine fill_merid_halo_cell_3d

  subroutine fill_zonal_halo_lon_2d(field, all_halo, left_halo, right_halo)

    real(r8), intent(inout) :: field(mesh%half_lon_lb:mesh%half_lon_ub,&
                                     mesh%full_lat_lb:mesh%full_lat_ub)
    logical, intent(in), optional :: all_halo
    logical, intent(in), optional :: left_halo
    logical, intent(in), optional :: right_halo

    integer i, j, m, n
    logical all_, left_, right_

    if (present(all_halo)) then
      all_ = all_halo
    else
      all_ = .true.
    end if
    if (present(left_halo)) then
      left_ = left_halo
    else
      left_ = .true.
    end if
    if (present(right_halo)) then
      right_ = right_halo
    else
      right_ = .true.
    end if

    if (all_ .or. left_) then
      ! left halo
      m = lbound(field, 1) - 1
      n = ubound(field, 1) - 2 * mesh%lon_halo_width
      do j = mesh%full_lat_ibeg+1, mesh%full_lat_iend-1
        do i = 1, mesh%lon_halo_width
          field(m+i,j) = field(n+i,j)
        end do
      end do
    end if

    if (all_ .or. right_) then
      ! right halo
      m = ubound(field, 1) - mesh%lon_halo_width
      n = lbound(field, 1) + mesh%lon_halo_width - 1
      do j = mesh%full_lat_ibeg+1, mesh%full_lat_iend-1
        do i = 1, mesh%lon_halo_width
          field(m+i,j) = field(n+i,j)
        end do
      end do
    end if 

  end subroutine fill_zonal_halo_lon_2d

  subroutine fill_zonal_halo_lon_3d(field, all_halo, left_halo, right_halo)

    real(r8), intent(inout) :: field(mesh%half_lon_lb:mesh%half_lon_ub, &
                                     mesh%full_lat_lb:mesh%full_lat_ub, &
                                     mesh%full_lev_lb:mesh%full_lev_ub)
    logical, intent(in), optional :: all_halo
    logical, intent(in), optional :: left_halo
    logical, intent(in), optional :: right_halo

    integer i, j, k, m, n
    logical all_, left_, right_

    if (present(all_halo)) then
      all_ = all_halo
    else
      all_ = .false.
    end if
    if (present(left_halo)) then
      left_ = left_halo
    else
      left_ = .false.
    end if
    if (present(right_halo)) then
      right_ = right_halo
    else
      right_ = .false.
    end if

    if (all_ .or. left_) then
      ! left halo
      m = lbound(field, 1) - 1
      n = ubound(field, 1) - 2 * mesh%lon_halo_width
      do k = mesh%full_lev_lb, mesh%full_lev_ub
        do j = mesh%full_lat_ibeg+1, mesh%full_lat_iend-1
          do i = 1, mesh%lon_halo_width
            field(m+i,j,k) = field(n+i,j,k)
          end do
        end do
      end do
    end if

    if (all_ .or. right_) then
      ! right halo
      m = ubound(field, 1) - mesh%lon_halo_width
      n = lbound(field, 1) + mesh%lon_halo_width - 1
      do k = mesh%full_lev_lb, mesh%full_lev_ub
        do j = mesh%full_lat_ibeg+1, mesh%full_lat_iend-1
          do i = 1, mesh%lon_halo_width
            field(m+i,j,k) = field(n+i,j,k)
          end do
        end do
      end do
    end if 

  end subroutine fill_zonal_halo_lon_3d

end module parallel_mod